//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : solvepsi_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef SOLVEPSI_KERNELS_H
#define SOLVEPSI_KERNELS_H

#ifdef USEDOUBLE
#define ZERO 0.0
#define QUART 0.25
#else
#define ZERO 0.0f
#define QUART 0.25f
#endif

namespace wavescomputer
{

namespace solvepsi
{

template <class T>
__global__ void kernel_phase1(T* bpsi,
                              const T* h0,
                              const T* D0,
                              const T* phi,
                              const T* SRCp,
                              const T* SppC, // <-- what is this?
                              const int ld,
                              const int Nx1,
                              const int Nx2,
                              const int dx1,
                              const int dx2,
                              const T pslow)
{
    int u = BORDER_WIDTH + bx * Bx + tx;
    int v = BORDER_WIDTH + by * By + ty;
    int tid = v * ld + u;  

    const T mult1 = (QUART / (dx1 * dx1));
    const T mult2 = (QUART / (dx2 * dx2));
    
    const T val_h0 = h0[tid];
    const T val_D0 = D0[tid];

    if (v - BORDER_WIDTH < Nx1  && u - BORDER_WIDTH < Nx2)
    {
        if (SppC[tid] > ZERO)
        {
            const T SpfE = (val_h0 + h0[tid + ld]) 
                         * (val_D0 + D0[tid + ld]) * mult1;
            const T SpfW = (val_h0 + h0[tid - ld]) 
                         * (val_D0 + D0[tid - ld]) * mult1;
            const T SpfN = (val_h0 + h0[tid + 1 ]) 
                         * (val_D0 + D0[tid + 1 ]) * mult2;
            const T SpfS = (val_h0 + h0[tid - 1 ]) 
                         * (val_D0 + D0[tid - 1 ]) * mult2;
            const T SpfC = -(SpfE + SpfW + SpfN + SpfS);

            bpsi[tid] = - SpfE * phi[tid + ld]
                        - SpfW * phi[tid - ld]
                        - SpfN * phi[tid + 1 ]
                        - SpfS * phi[tid - 1 ]
                        - SpfC * phi[tid]
                        - SRCp[tid] * pslow;
        }
        else
        {
            bpsi[tid] = - SRCp[tid] * pslow;
        }
    }
}

namespace reflectiveneumann
{

// Reflective Neumann boundary (North version)
template <class T>
__global__ void kernel_sp_rn_north(T* bpsi, 
                                   const T* stoN,
                                   const T* psi_in_N_ext, 
                                   const T* psi_in_N_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2 - 1); // not coalesced

    if (tid < Nx1) 
    {
        bpsi[tid2D] -= stoN[tid1D] * (psi_in_N_ext[tid1D] - psi_in_N_cen[tid1D]);
    }
}

// Reflective Neumann boundary (East version)
template <class T>
__global__ void kernel_sp_rn_east(T* bpsi, 
                                  const T* stoE,
                                  const T* psi_in_E_ext, 
                                  const T* psi_in_E_cen, 
                                  const int ld,
                                  const int Nx1,
                                  const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + Nx1 - 1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
        bpsi[tid2D] -= stoE[tid1D] * (psi_in_E_ext[tid1D] - psi_in_E_cen[tid1D]);
    }
}

// Reflective Neumann boundary (South version)
template <class T>
__global__ void kernel_sp_rn_south(T* bpsi, 
                                   const T* stoS,
                                   const T* psi_in_S_ext, 
                                   const T* psi_in_S_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH;             // not coalesced

    if (tid < Nx1) 
    {
        bpsi[tid2D] -= stoS[tid1D] * (psi_in_S_ext[tid1D] - psi_in_S_cen[tid1D]);
    }
}

// Reflective Neumann boundary (West version)
template <class T>
__global__ void kernel_sp_rn_west(T* bpsi, 
                                  const T* stoW,
                                  const T* psi_in_W_ext, 
                                  const T* psi_in_W_cen, 
                                  const int ld,
                                  const int Nx1,
                                  const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = BORDER_WIDTH * ld + BORDER_WIDTH + tid;               // coalesced

    if (tid < Nx2) 
    {
        bpsi[tid2D] -= stoW[tid1D] * (psi_in_W_ext[tid1D] - psi_in_W_cen[tid1D]);
    }
}

} // namespace reflectiveneumann

namespace sommerfeld
{

// Sommerfeld boundary (North version)
template <class T>
__global__ void kernel_sp_som_north(T* bpsi,
                                    const T* stoN,
                                    const T* psiN,
                                    const T* psi_in_N_ext,
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2 - 1); // not coalesced

    if (tid < Nx1) 
    {
        bpsi[tid2D] -= stoN[tid1D] * (psiN[tid1D] + psi_in_N_ext[tid1D]);
    }
}

// Sommerfeld boundary (East version)
template <class T>
__global__ void kernel_sp_som_east(T* bpsi,
                                   const T* stoE,
                                   const T* psiE,
                                   const T* psi_in_E_ext,
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + Nx1 - 1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
        bpsi[tid2D] -= stoE[tid1D] * (psiE[tid1D] + psi_in_E_ext[tid1D]);
    }
}

// Sommerfeld boundary (South version)
template <class T>
__global__ void kernel_sp_som_south(T* bpsi,
                                    const T* stoS,
                                    const T* psiS,
                                    const T* psi_in_S_ext,
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH;             // not coalesced

    if (tid < Nx1) 
    {
        bpsi[tid2D] -= stoS[tid1D] * (psiS[tid1D] + psi_in_S_ext[tid1D]);
    }
}

// Sommerfeld boundary (West version)
template <class T>
__global__ void kernel_sp_som_west(T* bpsi,
                                   const T* stoW,
                                   const T* psiW,
                                   const T* psi_in_W_ext,
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                                  // coalesced
    int tid2D = BORDER_WIDTH * ld + BORDER_WIDTH + tid;               // coalesced

    if (tid < Nx2) 
    {
        bpsi[tid2D] -= stoW[tid1D] * (psiW[tid1D] + psi_in_W_ext[tid1D]);
    }
}

} // namespace sommerfeld

} // namespace solvepsi

} // namespace wavescomputer

#endif // SOLVEPSI_KERNELS_H

